↓↓ Accès au projet ↓↓
Lien GitHub du projet

Introduction

Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des opérations statistiques avancées grâce au langage R. Ainsi, notre choix s’est porté sur des données statistiques relatives au prestige (score attribuée) de différentes professions canadiennes.

Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer une modélisation statistique sur notre jeu de données et produire un rapport.


Outils et environnement de travail

R et RStudio

R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.

La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.


R markdown

L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.


Outils complémentaires

  • tidyverse : Il s’agit d’une collection d’extensions relatives à la science des données, et permettant la manipulation des tableaux de données, l’import/export de données, la manipulation de variables ou, entre autres, la visualisation de données ;
  • data.table : Extension et manipulation avancée des tableaux de données ;
  • plotly : Comme ggplot2 (compris dans le tidyverse), ce package permet la visualisation de données, à la différence près qu’il renden ces graphiques interactifs ;
  • car : Appliquer une régression linéaire et des tests statistiques ;
  • hrbrthemes : Thèmes pour ggplot2.


Données

Notre jeu de données comporte 102 individus décrits par 6 variables :


head(data) %>% paged_table()
viz_edu <- data %>% 
  ggplot(aes(x = prestige, y = education, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(6, 16, 2)
  )

viz_edu <- ggplotly(viz_edu)

viz_inc <- data %>% 
  ggplot(aes(x = prestige, y = income, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(0, 20000, 10000)
  )

viz_inc <- ggplotly(viz_inc)

viz_women <- data %>% 
  ggplot(aes(x = prestige, y = women, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
      breaks = seq(25, 75, 25)
    ) +
    scale_y_continuous(
      breaks = seq(0, 100, 25)
    )

viz_women <- ggplotly(viz_women)

viz_census <- data %>% 
  ggplot(aes(x = prestige, y = census, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(2500, 7500, 2500)
  )
  
viz_census <- ggplotly(viz_census)

subplot(viz_edu, viz_inc, viz_women, viz_census, nrows = 2, margin = 0.07)

Régression linéaire

En statistiques et en apprentissage automatique, un modèle de régression linéaire est un modèle de régression qui cherche à établir une relation linéaire entre une variable, dite expliquée, et une ou plusieurs variables, dites explicatives.

Dans un premier temps, nous pouvons étudier rapidement d’éventuelles corrélations entre deux variables à l’aide d’un nuage de points. Ici, nous utilisons la fonction ggpairs afin de pouvoir en observer plusieurs en même temps.

generate_scatterM <- function(data, mapping){
  data %>% ggplot(mapping = mapping) + geom_point(size = 0.8) + 
    geom_smooth(method = loess, color = "red", size = 0.85, se = FALSE) +
    geom_smooth(method = lm, color = "blue", size = 0.85, se = FALSE)
}

scatterM <- data %>% ggpairs(columns = 1:5, lower = list(continuous = generate_scatterM))

scatterM <- ggplotly(scatterM)
scatterM

Avec les observations précédentes, nous pouvons déjà écarter certaines linéarités et choisir lesquelles sont intéressantes à étudier. Ici, nous reprenons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable education (qui reflète le niveau d’étude). Cette fois nous utilisons ggplot pour une meilleure représentation.

education <- data %>% 
  ggplot(aes(x = education, y = prestige)) +
  geom_point() +
  geom_smooth(method = loess, se = T) +
  geom_smooth(method = lm, color = "red", se = F) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(6, 16, 2)
  ) +
  scale_y_continuous(
    breaks = seq(20, 80, 20)
  )

scatter_edu <- ggplotly(education)
scatter_edu

La droite bleue est définie par la méthode des moindres carrés (MSE), il s’agit d’une droite de régression linéaire entre les variables education et prestige.

La courbe rouge indique la tendance globale entre ces deux variables, il s’agit d’une courbe de régression de type lowess. Les deux courbes extérieures représentent l’intervalle de confiance de cette courbe de régression.

On constate que que la droite de régression est presque toujours comprise dans l’intervalle de confiance, l’hypothèse de linéarite entre les variables education et prestige est donc acceptable.


Nous testons aussi la relation entre les variables income et prestige.

income <- data %>% 
  ggplot(aes(x = income, y = prestige)) +
  geom_point() +
  geom_smooth(method = loess, se = T) +
  geom_smooth(method = lm, color = "red", se = F) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(0, 25000, 5000)
  ) +
  scale_y_continuous(
    breaks = seq(20, 100, 20)
  )

scatter_inc <- ggplotly(income)
scatter_inc

Ici, en regardant la forme du lien entre la variable entre les deux variables, on s’aperçoit que la droite de régression suis bien moins l’intervalle de confiance de la courbe lowess. l’hypothèse de linéarité est alors plus critiquable.

Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction “lm”.

lin_model_edu <- lm(prestige ~ education, data = data)
lin_model_edu
## 
## Call:
## lm(formula = prestige ~ education, data = data)
## 
## Coefficients:
## (Intercept)    education  
##     -10.732        5.361

« Intercept » correspond ici à l’ordonnée à l’origine, le « b » de notre droite, et le « x » est la pente de la droite, ce qui correspond au « b » dans notre notation.

L’équation, de notre droite est donc \(y = -10.732 + 5.361x\).

bacf <- acf(residuals(lin_model_edu), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(lin_model_edu)))

lag_plot <- bacfdf %>%
  ggplot(aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(aes(xend = lag, yend = 0)) +
  geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +     
  geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
  theme_bw()

lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plot
durbinWatsonTest(lin_model_edu)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2752512       1.43917   0.002
##  Alternative hypothesis: rho != 0
qq_edu <- lin_model_edu %>%
  ggplot(aes(sample = rstandard(.))) +
  stat_qq_line(color = "red", linetype = "dashed") +
  stat_qq(size = 1.2) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(-2, 2, 1)
  ) +
  scale_y_continuous(
    breaks = seq(-3, 3, 1)
  ) +
  labs(title = "Normal Q-Q",
       x = "Theoretical Quantiles\nlm(prestige ~ education)", 
       y = "Standardized residuals")

qq_edu <- ggplotly(qq_edu)
qq_edu
shapiro.test(residuals(lin_model_edu))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lin_model_edu)
## W = 0.98065, p-value = 0.1406
lin_model_inc <- lm(prestige~income, data = data)
lin_model_inc
## 
## Call:
## lm(formula = prestige ~ income, data = data)
## 
## Coefficients:
## (Intercept)       income  
##   27.141176     0.002897
qq_inc <- lin_model_inc %>%
  ggplot(aes(sample = rstandard(.)))+
  stat_qq_line(color = "red", linetype = "dashed") +
  stat_qq(size = 1.2) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(-2, 2, 1)
  ) +
  scale_y_continuous(
    breaks = seq(-3, 3, 1)
  ) +
  labs(title = "Normal Q-Q",
       x = "Theoretical Quantiles\nlm(prestige ~ income)", 
       y = "Standardized residuals")

qq_inc <- ggplotly(qq_inc)
qq_inc
shapiro.test(residuals(lin_model_inc))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lin_model_inc)
## W = 0.97169, p-value = 0.02729
sc_loc_edu <- lin_model_edu %>%
  ggplot(aes(fitted(.), sqrt(abs(rstandard(.))))) + 
  geom_point(size = 1.2) +
  geom_smooth(method = loess, se = FALSE) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(30, 70, 10)
  ) +
  scale_y_continuous(
    limits = c(0, 1.5),
    breaks = seq(0, 1.5, 0.5)
  ) +
  labs(title = "Scale-location",
       x = "Fitted values\nlm(prestige ~ education)", 
       y = "sqrt |Standardized residuals|")

sc_loc_edu <- ggplotly(sc_loc_edu)
sc_loc_edu
ncvTest(lin_model_edu)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.6327545, Df = 1, p = 0.42635
res_v_fit_edu <- lin_model_edu %>%
  ggplot(aes(fitted(.), residuals(.))) + 
  geom_point() +
  stat_smooth(method="loess", se = FALSE) + 
  geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(30, 70, 10)
  ) +
  scale_y_continuous(
    limits = c(-30, 20),
    breaks = seq(-30, 20, 10)
  ) +
  labs(title = "Residuals vs. Fitted values",
       x = "Fitted values\nlm(prestige ~ education)", 
       y = "Residuals")
    
res_v_fit_edu <- ggplotly(res_v_fit_edu)
res_v_fit_edu
summary(lin_model_edu)
## 
## Call:
## lm(formula = prestige ~ education, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.0397  -6.5228   0.6611   6.7430  18.1636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -10.732      3.677  -2.919  0.00434 ** 
## education      5.361      0.332  16.148  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared:  0.7228, Adjusted R-squared:   0.72 
## F-statistic: 260.8 on 1 and 100 DF,  p-value: < 2.2e-16
confint(lin_model_edu) %>% kable()
2.5 % 97.5 %
(Intercept) -18.027220 -3.436744
education 4.702223 6.019533 
predictions <- tibble(education=c(9.21, 11.07, 14.62))
predict(lin_model_edu, newdata = predictions, interval = "confidence") %>% kable()
fit lwr upr
38.64170 36.58966 40.69374
48.61293 46.81134 50.41452
67.64405 64.52387 70.76423 
data$residuals <- residuals(lin_model_edu)
head(data) %>% kable()
education income women prestige census type residuals
13.11 12351 11.16 68.8 1113 prof 9.250875
12.26 25879 4.02 69.1 1130 prof 14.107621
12.77 9271 15.70 63.4 1171 prof 5.673573
11.42 8865 9.11 56.8 1175 prof 6.310758
14.62 8403 11.68 73.5 2111 prof 5.855950
15.64 11030 5.13 77.6 2113 prof 4.487854 
data$fitted <- fitted(lin_model_edu)
head(data) %>% kable()
education income women prestige census type residuals fitted
13.11 12351 11.16 68.8 1113 prof 9.250875 59.54913
12.26 25879 4.02 69.1 1130 prof 14.107621 54.99238
12.77 9271 15.70 63.4 1171 prof 5.673573 57.72643
11.42 8865 9.11 56.8 1175 prof 6.310758 50.48924
14.62 8403 11.68 73.5 2111 prof 5.855950 67.64405
15.64 11030 5.13 77.6 2113 prof 4.487854 73.11215 
vcov(lin_model_edu) %>% kable()
(Intercept) education
(Intercept) 13.520978 -1.1835055
education -1.183506 0.1102162 
pred_interval <- predict(lin_model_edu, interval = "prediction")
data <- cbind(data, pred_interval)
linear_regression <- data %>% 
  ggplot(aes(y = prestige, x = education)) +
  geom_point(size = 1.2) +
  geom_smooth(color = "red", method = "lm", fill = "red") +
  geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
  geom_line(aes(y =  upr), color = "red", linetype = "dashed") +    
  annotate("text", x = 9, y = 80, 
           label = paste0("prestige = ", coef(lin_model_edu)["(Intercept)"] %>% round(3), 
                          " + ", coef(lin_model_edu)["education"] %>% round(3), " * education")) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(6, 16, 2)
  ) +
  scale_y_continuous(
    breaks = seq(25, 75, 25)
  ) +
  labs(title = "Linear Regression",
       x = "Education", y = "Prestige")
  

linear_regression <- ggplotly(linear_regression)
linear_regression